W rozdziale tym przyjrzymy się jakie macierze generują opisane powyżej funkcje bazowe B-spline, i oszacujemy koszt faktoryzacji wynikającego z nich układu równań. W ten sposób możemy wskazać funkcje bazowe, które dają szybsze obliczenia za pomocą izogeometrycznej metody elementów skończonych. Ustalmy dla uproszczenia że rozwiązujemy problem projekcji jednowymiarowej. Szczegóły na temat obliczeń dwuwymiarowej projekcji zostały opisane w rozdziale pierwszym. Mamy wówczas następujące macierze zwane macierzami masowymi.
Załóżmy, że chcemy rozwiązać jednowymiarowy problem projekcji.
\( u(x) \approx BITMAP(x) \)
Szczegóły na temat obliczeń dwuwymiarowej projekcji zostały opisane w rozdziale pierwszym. Tutaj, dla uproszczenia załóżmy, że rozwiązujemy jednowymiarowy problem projekcji, nasza bitmapa (którą aproksymujemy w rozdziale pierwszym) jest teraz jednowymiarową płaską bitmapą.
Nasza jednowymiarowa bitmapa jest aproksymowana przez kombinację liniową funkcji B-spline
\( u = \sum_{i=1,...,N_x} u_{i} B^x_{i;p }(x) \)
Załóżmy, że mamy 16 elementów skończonych (16 przedziałów). Oczywiście nie oznacza to że \( N_x=16 \), ponieważ \( N_x \) to liczba funkcji bazowych. To, ile będzie funkcji bazowych, zależy od tego jak wyglądać będzie nasz wektor węzłów, innymi słowy od tego jakiego stopnia będą nasze funkcje bazowe, oraz czy i w których miejscach powtórzymy węzły pomiędzy elementami.
W celu wygenerowania układu równań z naszą macierzą do problemu projekcji, przemnażamy nasze równanie przez poszczególne B-spliny, które służą do uśredniania naszej relacji \( u(x)=BITMAP(x) \) zgodnie z rozkładem opisanym przez B-spline'y, i całkujemy żeby policzyć średnią
\( \int{\color{red}{u(x)}}B^x_{1;p}(x)dx=\int BITMAP(x)B^x_{1;p }(x)dx \\ \int{\color{red}{u(x)}}B^x_{2;p}(x)dx=\int BITMAP(x)B^x_{2;p}(x)dx \\ \vdots \\ \int{\color{red}{u(x)}}B^x_{N_x-1;p}(x)dx=\int BITMAP(x)B^x_{N_x-1;p}(x)dx \\\int{\color{red}{u(x)}}B^x_{N_x;p}(x)dx=\int BITMAP(x)B^x_{N_x;p }(x)dx \)
Wstawiamy naszą aproksymację w miejsce \( {\color{red}{u(x)}}=\sum_{i=1,...,N_x} u_{i} B^x_{i;p }(x) \) i dostajemy układ równań:
\( \int{\color{red}{\sum_{i=1,...,N_x} u_{i} B^x_{i;p}(x)} }B^x_{1;p}(x)dx=\int BITMAP(x)B^x_{1;p}(x) dx \\ \int{\color{red}{\sum_{i=1,...,N_x} u_{i} B^x_{i;p}(x)} }B^x_{2;p}(x)dx=\int BITMAP(x)B^x_{2;p}(x) dx \\ \vdots \\ \int{\color{red}{\sum_{i=1,...,N_x} u_{i} B^x_{i;p}(x)} }B^x_{N_x-1;p}(x)dx=\int BITMAP(x)B^x_{N_x-1;p}(x) dx \\ \int{\color{red}{\sum_{i=1,...,N_x} u_{i} B^x_{i;p}(x)} }B^x_{N_x;p}(x)dx=\int BITMAP(x)B^x_{N_x;p}(x) dx \)
Wyciągamy sumę przed całkę
\( \sum_{i=1,...,N_x}\int{\color{red}{ u_{i} B^x_{i;p}(x)} }B^x_{1;p}(x)dx=\int BITMAP(x)B^x_{1;p}(x) dx \\ \sum_{i=1,...,N_x} \int{\color{red}{u_{i} B^x_{i;p}(x)} }B^x_{2;p}(x)dx=\int BITMAP(x)B^x_{2;p}(x) dx \\ \vdots \\ \sum_{i=1,...,N_x} \int{\color{red}{u_{i} B^x_{i;p}(x)} }B^x_{N_x-1;p}(x)dx=\int BITMAP(x)B^x_{N_x-1;p}(x) dx \\ \sum_{i=1,...,N_x} \int{\color{red}{ u_{i} B^x_{i;p}(x) } }B^x_{N_x;p}(x)dx=\int BITMAP(x)B^x_{N_x;p}(x) dx \)
Możemy zapisać nasz układ równań w następującej postaci
\( \begin{bmatrix}\int B^x_{1;p}(x)B^x_{1;p}(x)dx & \int B^x_{1;p}(x)B^x_{2;p}(x)dx & \cdots & \int B^x_{1;p}(x)B^x_{N_x;p}(x)dx \\ \int B^x_{2;p}(x)B^x_{2;p}(x)dx & \int B^x_{2;p}(x)B^x_{2;p}(x)dx & \cdots & \int B^x_{2;p}(x)B^x_{N_x;p}(x)dx \\ \vdots & \vdots & \vdots & \vdots \\ \int B^x_{N_x-1;p}(x)B^x_{1;p}(x)dx & \int B^x_{N_x-1;p}(x)B^x_{2;p}(x)dx & \cdots & \int B^x_{N_x-1;p}(x)B^x_{N_x;p}(x)dx \\ \int B^x_{N_x;p}(x)B^x_{N_x;p}(x)dx & \int B^x_{N_x;p}(x)B^x_{2;p}(x)dx & \cdots & \int B^x_{N_x;p}(x)B^x_{N_x;p}(x)dx \\ \end{bmatrix} \begin{bmatrix} u_{1} \\ u_{2} \\ \vdots \\ u_{N_x-1}\\ u_{N_x}\\ \end{bmatrix} = \)
\( \begin{bmatrix} \int BITMAP(x)B^x_1(x)dx \\ \int BITMAP(x)B^x_2(x)dx \\ \vdots \\ \int BITMAP(x)B^x_{N_x}-1(x)dx \\ \int BITMAP(x)B^x_{N_x}(x)dx \\ \end{bmatrix} \)

Zauważmy, że układ równań jest rzadki ponieważ całki z iloczynów jednowymiarowych funkcji B-spline są niezerowe jedynie wtedy gdy funkcje bazowe występujące w całce nachodzą na siebie.

Dlaczego tak jest?

Całka to suma próbek wartości funkcji, które całkujemy. Nasza próbka to iloczyn dwóch funkcji \( i \)-tej oraz \( j \)-tej, \( B^x_{i;p}(x)B^x_{j;p}(x) \). Jeśli więc dla danego punktu w którym próbkujemy jedna z funkcji, na przykład \( i \)
-ta jest równa zero, mamy więc \( 0*B^x_{j;p}(x)=0 \). Podobnie jeśli druga funkcja \( j \)-ta jest równa zero w miejscu próbkowania, mamy wówczas \( B^x_{i;p}(x)*0=0 \).

Załóżmy teraz, że na naszych szesnastu elementach rozpinamy wektor węzłów
[0 0 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 16 16 16].
Uzyskujemy w ten sposób funkcje B-spline trzeciego stopnia, stosowane w izogeometrycznej metodzie elementów skończonych. Są one przedstawione na Rys. 2. Rysunek ten ilustruje również fakt, iż kolejne wiersze i kolumny w macierzy związane są z kolejnymi funkcjami B-spline, określonymi w wektorze węzłów. Zauważmy, że każda jednowymiarowa funkcja B-spline stopnia \( p \) nachodzi na \( 2p+1 \) innych jednowymiarowych funkcji B-spline. Na przykład B-spline'y stopnia trzeciego \( p=3 \) nachodzą na \( 2*3+1=7 \) sąsiednich B-spline'ów.
Niezerowe wartości w macierzy występują w tych wierszach i kolumnach, dla których odpowiadające im funkcje bazowe "nachodzą" na siebie (mają nie zerowe przecięcie dziedzin).
Nasz układ równań jest więc rzadki, oraz zawiera on siedem przekątnych. W szczególności dla naszego przypadku szesnastu elementów skończonych mamy \( N=16+3=19 \) funkcji bazowych B-spline trzeciego stopnia (w ogólności jest ich \( N_e+p \)). Nasza macierz ma więc rozmiar \( 19 \times 19 \) i zawiera 7 gęstych przekątnych. Jest to zilustrowane na Rys. 2.

Oszacujmy teraz koszt faktoryzacji macierzy dla bazy B-spline'ów trzeciego stopnia używanych w izogeometrycznej metodzie elementów skończonych, przedstawionej na Rys. 2. Uruchamiamy algorytm eliminacji Gaussa, który omija wszystkie zera w macierzy. Jest on opisany w rozdziale czwartym.
Koszt ten równy jest kosztowi faktoryzacji pojedynczego wiersza na podmacierzy o rozmiarze \( 4\times 4 \) przemnożony przez wiele wierszy macierzy, czyli 4 (koszt skalowania wiersza o czterech kolumnach) plus 3*3*2 (koszt odejmowania pierwszego wiersza od pozostałych, wraz z mnożeniem przez wartość z przekątnej). Operację tę powtarzamy dla 16 bloków macierzy (odpowiadających 16 elementom), a dla ostatniego bloku dokańczamy faktoryzację co ma koszt 2*2*2 (koszt odejmowania drugiego wiersza od pozostałych, wraz z mnożeniem przez wartość z przekątnej) + 2*1 (koszt odejmowania trzeciego wiersza od czwartego, wraz z mnożeniem przez wartość z przekątnej) + 1 (skalowanie ostatniego wiersza).Koszt całkowity to 4*16+3*3*2*16+2*2*2+2*1+1=64+18*16+8+2+1=363.

Rozważmy teraz przypadek, w którym powtarzamy węzły pomiędzy poszczególnymi elementami skończonymi. W ogólności jeśli chcemy uzyskać funkcje bazowe równoważne wielomianom Lagrange'a używanym w klasycznej metodzie elementów skończonych, musimy powtórzyć węzły \( p-1 \) razy, gdzie \( p \) to stopień bazy B-spline'ów. Jeśli więc stosujemy B-spliny trzeciego stopnia, musimy powtórzyć węzły 2 razy:
[0 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9 10 10 10 11 11 11 12 12 12 13 13 13 14 14 14 15 15 15 16 16 16 16]

Jak wygląda macierz naszego układu równań? Każda funkcja bazowa trzeciego stopnia rozpościera się na czterech przedziałach (używa pięciu węzłów, począwszy od 0 0 0 0 1 do 15 16 16 16 16. Ile mamy takich podprzedziałów w naszym wektorze węzłów? Jest ich 16*3+1=49. Nasza macierz ma więc wymiar \( 49 \times 49 \), i jest przedstawiona na Rys. 3.
Różnica macierzy właśnie uzyskanej dla bazy równoważnej wielomianom Lagrange'a w stosunku do macierzy uzyskanej wcześniej dla funkcji B-spline jest następująca. Macierz dla funkcji B-spline jest mniejsza, ale posiada gęstsze przekątne. Macierz uzyskana dla bazy równoważnej bazie Lagrange'a jest większa, ale za to posiada rzadsze przekątne.

Oszacujmy teraz koszt faktoryzacji macierzy dla bazy równoważnej bazie Lagrange'a używanej w standardowej metodzie elementów skończonych, przedstawionej na Rys. 3. Koszt ten równy jest kosztowi faktoryzacji pojedynczej podmacierzy o rozmiarze \( 49 \times 49 \) przez pełny algorytm eliminacji Gaussa, czyli 4+3+2 (koszt skalowania pierwszego, drugiego, trzeciego i czwartego wiersza od przekątnej do końca wiersza) plus 3*3*2 (koszt odejmowania pierwszego wiersza od pozostałych, wraz z mnożeniem przez wartość z przekątnej) plus 2*2*2 (koszt odejmowania drugiego wiersza od pozostałych, wraz z mnożeniem przez wartość z przekątnej) + 1*2 (koszt odejmowania trzeciego wiersza od czwartego, wraz z mnożeniem przez wartość z przekątnej). Koszt całkowity dla jednego bloku to 4+3+2+3*3*2+2*2*2+2*1=9+18+8+2=36. Mamy 16 elementów, więc koszt całkowity to 36*16=576.

Zobaczmy teraz co stanie się jeśli pomieszamy wielomiany Lagrange'a z funkcjami B-spline. Wprowadźmy wektor węzłów
[0 0 0 0 1 2 3 4 4 4 5 6 7 8 8 8 9 10 11 12 12 12 13 14 15 16 16 16 16],
w którym powtórzyliśmy węzły na grupach co cztery elementy. Uzyskamy wtedy funkcje bazowe i macierz przedstawioną na Rys. 4. W tym przypadku mamy 25 funkcji bazowych. Nasza macierz ma więc rozmiar \( 25\times 25 \), i posiada gęste bloki na przekątnych o rozmiarze 7 na 7. Bloki te poprzedzielane są separatorami. Jaki jest koszt faktoryzacji w tym przypadku?
Koszt ten równy jest kosztowi faktoryzacji pojedynczego wiersza na podmacierzy o rozmiarze \( 4\times 4 \) przemnożony przez wiersze macierzy w pojedynczym bloku, czyli 4 (koszt skalowania wiersza o czterech kolumnach) plus 3*3*2 (koszt odejmowania pierwszego wiersza od pozostałych, wraz z mnożeniem przez wartość z przekątnej). Operacje tą powtarzamy dla 4 podmacierzy w bloku (odpowiadających 4 elementom), i całość powtarzamy 4 razy (ponieważ mamy 4 bloki w macierzy). Całkowity koszt wynosi więc (4+3*3*2*4)*4=76*4=304.

Podsumujmy. Klasyczna metoda elementów skończonych dla 16. elementów i bazy równoważnej bazie Lagrange'a trzeciego stopnia wymaga 576 operacji zmienno-przecinkowych żeby sfaktoryzować macierz (mnożeń i dodawań / odejmowań). Izogeometryczna metoda elementów skończonych dla 16 elementów bazy równoważnej bazie B-spline'ów trzeciego stopnia wymaga 363 operacji zmienno-przecinkowych. Hybryda wymaga 304 operacji zmienno-przecinkowych.

Widzimy więc że nawet w prostym jednowymiarowym przykładzie hybryda pomiędzy bazą funkcji B-spline oraz bazą równoważną bazie Lagrange'a wygrywa w sensie kosztu obliczeniowego z czystą bazą B-spline'ów oraz z czystą bazą równoważną bazie Lagrange'a. Różnica ta będzie dużo większa (do dwóch rzędów wielkości) dla problemów dwuwymiarowych oraz trójwymiarowych.

Obliczenia za pomocą metody elementów skończonych, w których używa się baz mieszanych, w których stosuje się powtarzanie węzłów w wektorze węzłów nazywane są ulepszoną analizą izogeometryczną (ang. refined isogeometric analysis) [1]
Zauważmy również, iż wsadzając \( C^0 \) separatory (redukując ciągłość na granicach elementów) niejako zwiększamy moc aproksymacyjną używanej bazy funkcji. Tak więc dodawanie nowych funkcji i redukowanie ciągłości zwiększa dokładność aproksymacji.